% source code from https://zhuanlan.zhihu.com/p/590467893
% [t,y] = ode45(odefun,tspan,y0)
% [t,y] = ode45(odefun,tspan,y0,options)
% [t,y,te,ye,ie] = ode45(odefun,tspan,y0,options)

%%定义初始化条件，参数请参考原文：
x0 =[10*cos(2*pi/30); 10*sin(2*pi/30)];           % 初始状态
setGlobalx(x0);                                   % 设置初始X0
fopen('ETC.txt','w');                             % 打开文件ETC.txt，并清零，避免上次操作的影响 这两个文件自行新建在同一个M文件夹下面
fopen('norme.txt','w');                           % 打开文件norme.txt，并清零，避免上次操作的影响
%%使用Matlab ODE45对以上定义的两个微分方程进行求解：
dt=0.001
% [t1,x] = ode45(@ETC2007,[0:dt:15],x0);                % ode45 事件触发控制仿真 区间xspan：0-15 初始值：X0 
% [t2,y] = ode45(@system2007,[0:dt:15],x0);             % ode45 原系统仿真       区间xspan：0-15 初始值：X0

%%绘制结果部分：
[t1,x] = ode45(@ETC2007,[0:dt:15],x0);                % ode45 事件触发控制仿真 区间xspan：0-15 初始值：X0 
[t2,y] = ode45(@system2007,[0:dt:15],x0);             % ode45 原系统仿真       区间xspan：0-15 初始值：X0
figure(1);
plot(t2,y(:,1),'p green',t2,y(:,2),'* black');  % 两个系统对比 :t,x(:,1),'r',t,x(:,2),'c',t,
hold on
plot(t1,x(:,1),'red',t1,x(:,2),'blue');  
figure(2);
etc=dlmread('ETC.txt'); 
stem(etc);     
norme=importdata('norme.txt');
figure(3);
c = linspace(0,15,883); 
plot(c,norme);


function dx = ETC2007(t,x)
%%定义第二个微分方程：事件触发的微分方程：
xt = getGlobalx;
e = xt-x;
if norm(e) > 0.05 * norm(x) && t>0
   setGlobalx(x); 
   disp(t);
   dlmwrite('ETC.txt',1,'-append');
   dlmwrite('norme.txt',norm(e),'-append');

else
    dlmwrite('ETC.txt',0,'-append');
    dlmwrite('norme.txt',norm(e),'-append');

end
xt = getGlobalx;
A = [0 1; -2 3];  %系统参数矩阵
B = [0; 1];       %系统参数矩阵
K = [1 -4];       %反馈矩阵
u = K*xt;
dx=A*x+B*u;
end

function dy = system2007(t,y)
%%定义第一个微分方程：无事件触发的微分方程：
A = [0 1; -2 3];  %系统参数矩阵
B = [0; 1];       %系统参数矩阵
K = [1 -4];       %反馈矩阵
dy=A*y+B*K*y;
%u = k*x;
end


%%首先存储/读取上一次触发时刻状态变量的函数：

function setGlobalx(val)     %%%%存储
global x;
x = val;
end
function r = getGlobalx             %%%%%读取
global x;
r = x;
end